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We investigate the dynamical behavior of simple modules composed of two genes with two or three 
regulating connections. Continuous dynamics for mRNA and protein concentrations is compared to 
a Boolean model for gene activity. Using a generalized method, we study within a single framework 
different continuous models and different types of regulatory functions, and establish conditions 
under which the system can display stable oscillations. These conditions concern the time scales, 
the degree of cooperativity of the regulating interactions, and the signs of the interactions. Not 
all models that show oscillations under Boolean dynamics can have oscillations under continuous 
dynamics, and vice versa. 



I. INTRODUCTION 

Soon after the first repressor protein had been discov- 
ered by Francois Jacob and Jacques Lucien Monod in 
1961 [1], theoretical work on gene regulation started [2]. 
Typically, a realistic modelling of gene regulation systems 
includes rate equations for the concentrations of the par- 
ticipating macromolecules, i.e., mRNA and proteins. In 
vivo experiments permitted to obtain quantitative data 
for regulatory processes and their kinetic parameters, and 
they provided important insight into regulatory dynam- 
ics. Thus, Elowitz and Leibler [3] designed and con- 
structed a synthetic network out of three transcriptional 
repressor systems to build an oscillating network, termed 
the repressilator, in Escherichia coli. In recent years, the 
field of systems biology has emerged, which aims at a 
quantitative description of cell behavior and basic dy- 
namic processes, and which permits to analyze systems 
such as gene regulation networks. 

Besides detailed quantitative approaches, the general 
features of regulatory and signaling processes in living 
cells also gave rise to a minimalist dynamical description 
as Boolean networks, where the state of each gene is ei- 
ther "on" or "off" [4, 5]. Such a description is particularly 
useful when dealing with large networks [6] , because it re- 
duces the huge complexity of the continuous system with 
its many differential equations and parameters to a prob- 
lem of logical structure which is easier to understand. It 
permits to study generic features of entire classes of sys- 
tems [4], or to reproduce the correct sequence of events 
in gene regulation networks that must function reliably, 
such as cell cycle dynamics [7]. 

So far, little is known about the general conditions un- 
der which a Boolean simplification gives a realistic pic- 
ture of the dynamics in gene regulation networks. In 
contrast to Boolean networks, ordinary differential equa- 
tions (ODEs), which model the switch- like dynamics of 
genes by using sigmoidal Hill functions, can include more 
detailed information about transcription and translation 
processes to evaluate the time course of the gene expres- 
sion patterns. Depending on the parameter values, such 
models can show oscillating behavior or a stable fixed 
point. Even for small systems of only two genes, the 
seems to be no simple relation between Boolean and con- 



tinuous models. Widder et al. [8] and Polynikis et al. [9] 
studied a two-gene activator-inhibitor network and inves- 
tigated in detail the conditions under which a Hopf bifur- 
cation occurs, which leads to oscillating behavior. They 
found that oscillatory behavior is exhibited by the two- 
gene model only if the Hill coefficents are above a certain 
threshold, and that the system can be driven through the 
bifurcation by varying the time scales of the mRNA and 
proteins [9]. The Boolean version of this system always 
shows oscillatory behavior. 

Similar results can be found in Del Vecchio [10], who 
included an additional self-input of one gene and var- 
ied the time-scale difference between the activator and 
the respressor and time-scale difference between the pro- 
tein and mRNA dynamics by using bifurcation analysis. 
They also obtained richer dynamical behavior incorpo- 
rating mRNA dynamics and could define a parameter 
space that leads to stable limit cycles. 

In this paper, we present a general and comprehen- 
sive investigation of the two-gene network. Gross and 
Feudel [11] developed a method of generalized models, 
which allows to investigate the stability of fixed points 
and the occurrence of bifurcations in dependence of gen- 
eral features of the system, without the need to specify 
the steady state or the regulatory functions. This ap- 
proach enables us to unite all previous studies of this 
system within one framework, and to include also those 
situations that had not been studied before. This permits 
us to identify the main differences between the dynam- 
ical behavior and the attractor patterns of Boolean and 
continuous models. 

The paper is structured as follows. In Section II, we 
introduce the dynamical equations and the generalized 
method used for the analysis. Section III presents the 
conditions for the occurrence of oscillations in the con- 
tinuous model and compares these to the Boolean model. 
Finally, we discuss and compare our findings to previous 
studies in Section IV. 



II. MODEL 

Gene expression is the process by which genetic infor- 
mation is ultimately transformed into working proteins. 
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The main steps are transcription from DNA to mRNA, 
translation from mRNA to linear amino acid sequences 
and folding of these into functional proteins [12]. A cer- 
tain class of proteins, called transcription factors, can 
bind to the DNA to regulate the rate at which their tar- 
get genes are transcribed into mRNA. Gene regulation 
thus involves a network of macromolecules that mutually 
influence each other. The production of proteins and 
mRNA is balanced by degradation and dilution [13]. 



a b mRNA a >- protein a 

protein h ~ mRNA b 

FIG. 1: An example of a simple gene regulatory network con- 
sisting of two genes a and b regulating each other, and an 
additional self-input of gene a. The network consists of four 
molecular species: mRNA a and b, and proteins a and b. 
Dotted lines label translation processes, continuous lines tran- 
scription. Protein a regulates the production of mRNA b via 
activation or inhibition. Protein a and b together regulate the 
production of mRNA a via a two-dimensional input function. 



In this paper, we consider a simple gene regulatory 
network, consisting of two genes a and b regulating each 
other and an additional self-input of gene a. The concen- 
tration of mRNA produced by gene i is denoted by Ri, 
while the corresponding protein concentration is denoted 
by Pi, for i = a,b. Translation of mRNA and degradation 
of mRNA and proteins are supposed to occur with con- 
stant rate. Therefore, the ordinary differential equation 
system describing the reaction kinetics is [9] 

R a = m a Fi(P a ,P b ) - j a R a 
R b = m b F 2 (P a ) - ibRb 

P a = U a Ra - 5 a P a 

P b = uj b R b - 5 b P b 

F\ and F 2 depend on the concentrations of the regu- 
latory proteins. Regulation by only one protein (such as 
regulation of mRNA b by protein a) is usually modelled 
by a monotonically increasing sigmoidal-shaped function 
when the protein is an activator, and by a decreasing 
function when the protein is an inhibitor. Experimental 
evidence [14] suggests the usage of Hill functions [15]. 
Thus, the function F 2 (P a ) is either the Hill function for 
activation 

F + {P a ,k a ,n a ) = ° , (2) 
or the Hill function for inhibition, 

p« 

F-(P a ,k a ,n a ) = l-F+(P a ,k a ,n a ) = ° p » o . (3) 



Ri Concentration of transcribed mRNAs 

Pi Concentration of translated proteins 

rtii Maximal transcription rates 

u)i Translation rates 

ji mRNA degradation rates 

Si Protein degradation rates 

ki Expression threshholds 

m Hill coefficients 

TABLE I: Notation 

k a is the activation coefficient or expression thresh- 
old that defines the concentration of protein a needed 
to significantly activate expression. The parameter n a , 
called Hill coefficient, controls the steepness of the Hill 
function. The larger n a , the more step- like is the regula- 
tory function. Biologically, n a is related to the molecular 
binding mechanism: it is the number of proteins required 
for saturation of binding to the DNA [8] and is therefore 
also called cooperativity coefficient. The Hill function can 
be considered to be the probability that the promoter re- 
gion is bound, averaged over many binding and unbinding 
events of proteins i [16]. 

The combined effect of multiple transcription factors 
is described by using multi-dimensional input functions. 
An example for deriving such a function from gene ex- 
pression measurements was given by Setty et al. [17], 
who used the lacZYA operon of Escherichia coli. The 
function F 2 (P a ,P b ) can, for instance, integrate an acti- 
vator P a and a repressor P b [16]. If the activator and the 
inhibitor bind to the promotor independently, there are 
four binding-states of the promotor: unbound, bound to 
either P a or P b , or bound to both proteins. Transcription 
occurs mainly in the case that the activator binds the pro- 
moter and the repressor does not, resulting in a P a AND 
NOT P b input function. This is one of the 16 possi- 
ble Boolean functions for two inputs, and the simplest 
gene regulation models use such Boolean functions. The 
16 Boolean functions can be represented in P a -Pf,-space 
as so-called BooleCubes. Whereever we need to specify 
the functional form of F 2 (P a , P b ), we will use a continu- 
ous generalization of BooleCubes to so-called HillCubes, 
which is based on Hill functions and was suggested by 
Wittmann et al. [18]. Thus, the HillCube of the input- 
function F 2 {P a ,P b ) = P a AND NOT P b is given by the 
expression F 2 {P a ,P b ) = F + (P a ,k a ,n a ) ■ F~(P b ,k b ,n b ) 
(shown in Fig. 2). The other HillCubes are constructed 
in analogous manner by taking sums over the appro- 
priate products of Hill functions. For instance, the 
XOR function is written as F 2 (P a , P b ) = F + (P a , k a ,n a ) ■ 
F~(P b ,k b , n b ) + F~(P a , k a ,n a ) ■ F+(P b , k b ,n b ). 

In order to investigate the conditions under which 
model (1) can display oscillations, we will perform a 
linear stability analysis of the fixed points, and deter- 
mine whether a Hopf bifurcation can occur. The fixed 
point condition is R a = R b = P a = Pb = 0, leading 
to R a * = ^P a * and R b * = ^P b *, and to the implicit 
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FIG. 2: HillCube of the input-function F = P a AND NOT P b . 
The concentration of transcribed mRNA is projected to the 
two-dimensional surface. In light-colored regions, the con- 
centration of activator, P a , is high, and the concentration of 
inhibitor, Pi, is low, so mRNA a is transcribed significantly. 



equations 



P* 



to a m a 
uj b m b 



(4) 



The Jacobian J at this fixed point can be written as 



Since we intend to investigate model (1) for various 
types of functions F 1 (P a ,Pi ) ) and of parameter values, 
we use the description in terms of generalized models, as 
introduced by Gross and Feudel [11]. This method al- 
lows to investigate the occurrence of different dynamical 
regions without the need to specify the precise form of 
interaction terms and equilibrium concentrations. This 
method has already been successfully applied to ecologi- 
cal food webs [11, 20, 21], socio-economic models [11] and 
metabolic networks [22, 23]. 

With the assumption that there exists at least one fixed 
point with positive concentrations R a * , Rb* , P a * \ Pb* , one 
can define normalized state variables n — -jp?,Pi = -p-V 

and normalized functions fj (pi ) = Fl ^p. ^ , with the ab- 
breviated notation F* = Fj(P*). The fixed point values 



then simply are r* 



Pa = Pb = 1 • Rewriting equa- 



tions (1) in terms of the normalized variables, we obtain 

-h(Pa,Pb) ~ l a r a 



P ! 
m b F£ ~ 



nr. 



h{pa) -ibn 



P",', 



Pa = pf W r a - S a p a 



Pi 



Pb 



UbTb ~ SbPb- 



(6) 



J = 




translation 



transcription 
protein degradation 



(5) 



Based on the eigenvalues of the Jacobian, one can de- 
termine parameter values at which bifurcations occur. 
The presence of a single zero eigenvalue indicates the 
emerge or destruction of fixed points or the exchange 
of stability properties between two fixed points. Bifur- 
cations of this type are the saddle node bifurcation, the 
transcritical bifurcation, and the pitchfork bifurcation. 
This type of bifurcations can be identified by the fact 
that the determinant of the Jacobian vanishes at the bi- 
furcation point. A second type of local codimension-1 
bifurcations is the Hopf bifurcation, characterized by the 
presence of two purely imaginary eigenvalues of the Ja- 
cobian. If the system is at a stable fixed point before the 
bifurcation, the steady-state becomes unstable at the bi- 
furcation, and if the bifurcation is supercricital, a limit 
cycle emerges in the vicinity of the unstable fixed point. 
Thus, a Hopf bifurcation implies a transition from sta- 
tionary to periodic motion. A method that yields analyt- 
ical test functions for the Hopf bifurcation is the method 
of resultants (described in detail in [19]). 



Introducing the four parameters 

a r = m a F*/R* a = j a , 
I3 r = m b FZ I PI = la , 
a p = R* a uj a /P* = 5 a , 
f3 p = Riuj b /P b = lb, 

we can rewrite our model as 

r a = OL r {fxip a ,Pb) - T a ) 

^b = Pr{h{Pa) - r b ) 
Pa = OL p {r a - pa) 

Pb = I3 p (r b - Pb). 



(7) 



(8) 



The Jacobian at a fixed point of this set of equations 



is 



J = 



«1, 




Throughout this paper, the terms with j £ {1,2} 
and i G {a, b} will be referred to as new variables fjPi, 
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which are called exponent parameters [11]. Positive val- 
ues of fjPi indicate that protein i is considered to be 
an activator. Negative values indicate inhibition, and if 
fjPi vanishes, protein i has no influence on the regulatory 
function Fj . Exponent parameters measure the degree of 
nonlinearity of a process in the steady-state [20]. We can 
interpret f2P a as an indicator of the level of the binding 
cooperativity of the protein, and therefore of the cooper- 
ativity coefficient n a . Using an activating Hill function, 
we obtain 



/(?«) 



F + (P a * Pa ) 

F+(P*) 



P*n a 



Pa 



p*n a n a , rn 
r a Pa i^a 



and 



df{p a ) 



d Pa 



Therefore, the exponent parameters fjPi are restricted 
to the interval [0,n], if protein i is an activator. In the 
same way, we find that fjPi G [— n, 0] if protein i is a 
repressor. 

Apart from the exponent parameters, the Jacobian also 
depends on the scale parameters [11] a r , a p , /3 r , and 
j3 p , which are a measure of the (inverse) time scales of 
the four concentrations. We assume identical ratios r = 
— = jt between the mRNA and protein dynamics for 
both genes. A large value of r means a quasi steady-state 
assumption for the mRNA, where the mRNA transients 
die out quickly. 

In the following, we will investigate the properties of 
the model in dependence of r and the three exponent 
parameters fj P i. 



III. RESULTS 

A. Model with only two connections 

We first consider the case that F± is independent of P a , 
leading to the following simplifed network 



The regulatory functions F\ and F^ can either have 
the same sign (activator-activator or inhibitor-inhibitor) 
or different signs (activator-inhibitor complex). In a 
Boolean model, an activating interaction is implemented 
by the Boolean function "copy", and a repressing inter- 
action by "invert". In the Boolean version, the model 
has only four possible states, 00, 01, 10, and 11. For 
the activator-activator or inhibitor-inhibitor case, the 
Boolean model has two fixed points and a cycle. The 
cycle alternates between 01 and 10 for the activator- 
activator, and between 00 and 11 for the inhibitor- 
inhibitor case. For the activator-inhibitor situation, the 



Boolean model has a cycle that involves all four states: 
00 -> 10 ->■ 11 -> 01 -> 00. We will see that the con- 
tinuous model can display a cycle only in the activator- 
inhibitor situation, and only if the parameters are in the 
appropriate range. 



-10 -5 

f2pa 
(a)r = 1 



-10 -5 5 

f2pa 
(b)r = 10 



f2pa 
(c)r = 50 



FIG. 3: Values of the exponent parameters fipt and fiVa £ 
{ — 10, 10}, for which the fixed point is unstable and has a com- 
plex conjugate pair of eigenvalues with a positive real part. 
The time scale ratio r between mRNA and protein dynamics 
increases from left to right. 

Figure 3 shows the regions in generalized parameter 
space where a fixed point is unstable against oscillations. 
The Hopf bifurcations, which occur at the boundary be- 
tween the white and the gray area, can be determined di- 
rectly from the characteristic polynomial of the Jacobian, 
which gives rise to only two eigenvalues. Oscillations can 
only occur when the exponent parameters fiPb and faPa 
have different signs, i.e., for the activator-inhibitor sys- 
tem. Furthermore, the product of the Hill coefficients 
n a and rib must be large enough. For larger time-scale 
ratio r between mRNA and protein dynamics, the expo- 
nent parameters must be larger to obtain a Hopf bifur- 
cation. For r = 50, where the mRNA dynamics is quasi 
always in a steady-state, the activator-inhibitor network 
shows no oscillatory dynamics for exponent parameters 
fiPb, jiPa £ {!)••) 4}, corresponding to cooperativity co- 
efficients ni £ {1,..,4}. Even though the generalized 
method uses normalised variables, the results can easily 
be compared to numerical simulations. Figure 4 shows 
mRNA and protein concentrations for different values of 
Ui (rows) with r = 1 (first column) and r — 50 (sec- 
ond column). As the separation of time scales between 
mRNA and proteins becomes larger, the oscillation fre- 
quency increases, but the amplitude of the oscillations 
decreases (n^ = 10) or even vanishes (rij < 10). 

In contrast to the Boolean model, where the state 
of a gene jumps instantaneously from "off" to "on" and 
back, the concentrations in the continuous model always 
change smoothly, even in the limit n — > oo, where the 
functions F\ and F2 become step functions. For this rea- 
son, the oscillation that occurs in the Boolean activator- 
activator (or inhibitor-inhibitor) system does not occur 
in the continuous model. Only the two fixed points that 
are also present in the Boolean model occur in the con- 
tinuous model. For the same reason, the oscillation of the 
Boolean activator-inhibitor system can occur in the con- 
tinuous model only when F\ and F2 are steep enough in 
order to drive the concentrations sufficiently fast through 
the intermediate values. Otherwise the system settles at 
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the (only) fixed point, which is found at intermediate 
concentration values. 



being stable. Figure 6 shows a cross section at f2p a 
for r = 1. 
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time scale ratio r = 1, Hill coefficient n = 3 
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time scale ratio r = 50, Hill coefficient n = 3 
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FIG. 5: Regions in parameter space where the fixed point has 
a complex pair of unstable eigenvalues, with the other two 
eigenvalues being stable. Time scale ratio r between mRNA 
and protein dynamics increases in the sequence of graphs. 
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FIG. 4: mRNA and protein concentrations of the network of 
one activator and one inhibitor for different values of m and 
r (equations (1) with Fi(P a ,P 6 ) = F+(P a ) and F 2 {P a ) = 
F~(P a )). The plots in the first column correspond to r = 1, 
the plots in the second column to r = 50, thus mRNA dynam- 
ics are 50 times faster, approximating the mRNA steady-state 
assumption. The different rows stand for different Hill coeffi- 
cients: rii = 2, 3, 5, 10. With incresing m, the amplitude of the 
oscillation increases and stabilizes. As the separation of time 
scales between mRNA and proteins becomes large (r = 50), 
the oscillation's frequency increases, but the amplitude of the 
oscillations gets smaller (rii = 10) or even vanishes. (Param- 
eter values: rrii = 2.0, ki — 0.5, 7* = 5, = u>i = 1.0) 



B. Model with three connections 

Figure 5 shows the generalization to the case f\p a ^ 
of Fig. 3. The Hopf bifurcations were calculated using 
the method of resultants [19]. Gray areas show regions in 
parameter space where the fixed point has a complex pair 
of unstable eigenvalues, with the other two eigenvalues 
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FIG. 6: Cross section through Fig. 5 (a) at f 2 p a = 3.0. 

One sees again that a faster mRNA dynamics leads to a 
smaller oscillatory region in parameter space. Just as for 
the activator-inhibitor network, and fcp a must have 
different signs to obtain a Hopf bifurcation, correspond- 
ing to opposed regulatory functions. The additional self- 
input of gene a can be activating or repressing for r = 1 , 
but must be activating for larger ratio of time scales (e.g. 
r = 50) to obtain oscillations. Larger values of the ex- 
ponent parameters, and therefore larger Hill coefficients 
m, can make up for a large ratio of time scales between 
mRNA and protein dynamics, as in the system without 
additional self input. However, biologically realistic val- 
ues of rii are in the range rii = 1 — 4 [16]. 

In order to compare the continuous model with the 
Boolean model, we now specify the functions f\ and to 
be HillCubes. Out of the 16 possible Boolean functions of 
two variables for F\(P a , Pb), only 10 actually depend on 
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the values of both variables. We discuss in the following 
these 10 functions and compare the Boolean dynamics 
with the dynamics due to HillCubes. We restrict our- 
selves to the case that a activates b, i.e., f2p a > 0. The 
case that a inhibits b can always be mapped on the first 
case in a Boolean model by inverting the states of nodes 
and changing the sign of the appropriate connections. 

The Boolean functions a AND b and a OR b give rise to 
the fixed points 00 and 11 in the Boolean model. Using 
HillCubes, we obtain fipt > for both functions, which 
means that the continuous model cannot have a Hopf 
bifurcation. Whether the continuous model has one or 
two stable fixed points, depends on the parameter values. 

The Boolean functions NOT a AND b and NOT a OR 6 
give rise to one fixed point (00 and 11, respectively) and 
one cycle involving the states 01 and 01 in the Boolean 
model. Using HillCubes, we obtain fipi, > for both 
functions, which means that the continuous model can- 
not have a Hopf bifurcation. This situation is analogous 
to the activator-activator loop, where the cycle present 
in the Boolean model does not occur in the continuous 
model. 

The Boolean functions a AND NOT b and a OR NOT b 
give rise to one fixed point (00 and 11, respectively), 
which is a global attractor, in the Boolean model. Us- 
ing HillCubles, we obtain fiph < and fip a > 0. 
The signs of the exponent parameters are such that a 
Hopf bifurcation is possible. When the concentration of 
transcribed mRNA is projected onto a two-dimensional 
surface (Fig. 2), one can see that the two functions 
a AND NOT b and a OR NOT b produce significant 
mRNA concentrations in the shaded area of Fig. 6. An 
optimal condition for a Hopf bifurcation are intermedi- 
ate values of the mRNA concentration, since the function 
F% is steep in this region. Figure 7 shows regions in pa- 
rameter space, where F2(P a ) is steep, for the function 
a AND NOT b (and for two other functions to be dis- 
cussed below). 

In order to determine the parameter regions that place 
the system in the shaded area of Fig. 6, we varied the pa- 
rameters rii, ki and r at fixed protein concentrations P* 
and -P h * . The values of P* and P h * were chosen such that 
the system is placed in the light-colored areas of Fig. 7, 
where oscillations are most likely. The result is shown 
in Fig. 8, and numerical simulations are shown in Fig. 9: 
There exists a complex pair of unstable eigenvalues for 
iii = 2 and r = 1. For ki = 0.35 we obtain sustained os- 
cillations. However, for this set of parameter values, the 
system is close to a bifurcation: at ki = 0.34 the periodic 
orbit grows until it collides with the fixed point 0. 

The Boolean functions a NOR b and a NAND b give 
rise to a cycle, which is a global attractor, in the Boolean 
model. Using HillCubles, we obtain fipi, < and 
fiPa < 0. The signs of the exponent parameters are 
such that a Hopf bifurcation is possible. The projections 
of the functions on a two-dimensional surface and the pa- 
rameter regions that place the system in the shaded area 
of Fig. 6 are again shown in Figs 7 and 8. The results are 




(a)P a ANDNOT P b {b)P a NOR P b 




0.0 0.2 0.4 0.6 0.8 1.0 



(c)P a XOR P b 

FIG. 7: Two-dimensional input functions F — 
P a AND NOT P b , F = P a NOR P b and F = P a XOR P b 
for ki — 0.5 and rii — 3, projected to the two-dimensional 
surface. The light-colored areas mark those regions in pa- 
rameter space where F2(P a ) is steep, implying intermediate 
values of the mRNA a concentration. 

compared to numerical simulations, shown in Fig. 10(b). 
In contrast to the previous case, the parameter region 
where oscillations are possible is much smaller. Oscilla- 
tions are only possible when mRNA dynamics is suffi- 
ciently slow and when the Hill coefficient n\, is above 2. 

Finally, we consider the Boolean functions a XOR b 
and a XNOR b. In contrast to the other functions dis- 
cussed so far, these functions are not canalyzing but 
change their output value whenever one of the input val- 
ues is changed. In the Boolean model, these functions 
give rise to one fixed point and to one cycle that com- 
prises the three remaining states. The signs of fip b and 
fiPa in the continuous model can now be positive or neg- 
ative, depending on the parameter values. An oscillation 
is thus possible in principle. Based on Fig. 7 we chose 
fixed point concentrations likely to induce oscillations to 
calculate the surface of the Hopf bifurcation. But in con- 
trast to the two sets of functions discussed in the previous 
section, we did not find an oscillatory region. 

IV. DISCUSSION AND CONCLUSIONS 

By using the method of generalized models to study the 
dynamics of simple two-gene regulatory network compo- 
nents, we could establish general conditions for the occur- 
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FIG. 9: mRNA and protein concentrations of the two-gene 
module with additional self-input and Fi = P a AND NOT P b 
(equations (1)) for rii — 2 and r = 1. For ki = 0.35 exists a 
complex pair of unstable eigenvalues, and we obtain sustained 
oscillations. At ki — 0.34, the oscillation has become unsta- 
ble, and the periodic orbit grows until it collides with the fixed 
point 0. (Parameter values: mi = 2.0, 7; = Si — u)i = 1.0) 
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time scale ratio r = 1, Hill coefficient n = 
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FIG. 8: Hopf bifurcation surface of the two-dimensional input 
functions F = a ANDNOT b and F = a NOR b. Within 
the surfaces, the fixed point has a complex pair of unstable 
eigenvalues, with the other two eigenvalues being stable. Us- 
ing generalized functions / of F = a ANDNOT 6 at a fixed 
protein concentration P£ = P b * = 0.5 and F — a NOR b 
at P* = 0.4, P b — 0.6, we varied m, ki and r. The choice 
of P* and P b is done with the help of Fig. 7, to be in an re- 
gion in parameter space which is likely to induce oscillations. 
The input function F — a XOR b does not provide a Hopf 
bifurcation. 



rence of Hopf bifurcations, which give rise to oscillatory 
dynamics. Apart from the signs of the regulatory inter- 
actions, the only relevant parameters in this general de- 
scription are the Hill coefficient and the time scale ratio 
between mRNA and protein dynamics. By comparing 
the different types of interactions to two-node Boolean 
models, we found that the occurrence of a cycle in the 
Boolean model is neither necessary nor sufficient for the 
occurrence of an oscillation in the continuous model. 

Our results combine and generalize the findings of 
several previous studies of such systems. The studies 
by Widder et al. [8] and Polynikis et al. [9] were fo- 
cused on a two-gene system with only two connections, 
and they found that oscillations can only occur in the 
activator-inhibitor case, and only for large Hill coeffi- 



time scale ratio r = 1, Hill coefficient n = 10 




} r = 10, Hill coefficient n = 10 

mRNA a 
mRNAb 
Protein a 

. , , , Protein b 



FIG. 10: mRNA and protein concentrations of the two-gene 
module with additional self-input and Pi = P a NOR P b 
(equations (1)) for different values of rn and r. The plots in 
the first column correspond to r — 1, the plots in the second 
column to r = 10, thus mRNA dynamics are 10 times faster, 
corresponding to Fig.5(b). The different rows stand for differ- 
ent Hill coefficients: rii — 2,3,5. With incresing m, the am- 
plitude of the oscillation increases and stabilizes. As the sep- 
aration of time scales between mRNA and proteins becomes 
larger (r = 10), the oscillations vanish for small values of m. 
(Parameter values: rrii = 2.0, ki = 0.5, 7i = 5i = u)i = 1.0) 



cients n > 3 [8], or more precisely for n a ■ rib > 4 [9]. 
(Choice of other parameters: 7; = Si = cj; = 1.0.) 
Polynikis et al. [9] found that the system can be driven 
through the Hopf bifurcation by varying the time scales 
of the mRNA and the proteins. Widder et al. [8] empha- 
size that the domain in parameter space that contains a 
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limit cycle becomes larger with increasing cooperativity 
as expressed by higher Hill coefficients. All these findings 
are contained concisely in our Fig. 3. 

The two-gene system with three connections was stud- 
ied by Del Vecchio [10], for the case fipt, < 0, f\p a > 0, 
HVa > 0, by using generalized Hill functions. The ex- 
tensive bifurcation analysis shows that oscillations occur 
over a larger parameter range when the mRNA dynamics 
is slower. This result is a special case of our very general 
result shown in Fig. 5. 

Norrell et al. [24], studied simple rings of four genes, 
and rings where one node has an additional self-input. 
They included a time delay to model translation, in- 
stead of including additional equations for the mRNA 
concentrations. They found that the time delay must be 
sufficiently large for oscillations to occur, which agrees 
with our finding that r must be sufficiently small. Fur- 
thermore, they found that continuous systems can ex- 
hibit stable oscillations in cases where Boolean reasoning 
would suggest otherwise, and that on a ring periodic cy- 
cles of Boolean systems do not exist in continuous sys- 
tems when the oscillation is not stable against fluctua- 
tions in the update time. This finding is a generalization 
of the result described above that the two-gene activator- 
activator system has no oscillations. In a different publi- 
cation [25], they point out the importance of the detailed 
form of the continuous functions for obtaining nontrivial 
dynamical patterns. This finding is related to our finding 
that the Hill coefficents rii must be sufficiently large for 
oscillations to occur. 

Mochizuki [26] investigated random networks with 
larger numbers of nodes. He found that the number of 



different steady states increases with the number of self- 
regulatory genes. He furthermore found that many of the 
periodic oscillations observed in the Boolean network are 
not present in the continuous model of gene regulatory 
networks and therefore the predictions of Boolean models 
can become unrealistic or too complex for larger networks 
when compared to those of the corresponding ODE mod- 
els. The latter finding is conform with the findings for 
the simple two-gene model. 

When all these investigations are taken together, there 
appears to be no simple criterion for deciding whether 
a periodic dynamical behavior in a Boolean model has 
an equivalent in a continuous model. In the examples 
investigated in this paper, there are two cases where the 
periodic dynamics of the Boolean model is "reliable" in 
the sense that fluctuations in the update times of the 
two nodes do not destroy the dynamical cycle. The first 
case is the activator-inhibitor system, the second is the 
NOR and NAND function. The corresponding contin- 
uous model has periodic oscillations whenenver the Hill 
coefficients are large enough and the time scale of the 
mRNA is slow enough. However, a reliable oscillation 
in the Boolean model is not a necessary requirement for 
an oscillation in the continuous model, as shown for the 
XOR and XNOR functions (where the oscillation in the 
Boolean model is not reliable), and for the AND NOT 
and OR NOT functions, where the Boolean model has a 
global fixed point. 

Certainly, more research is needed to establish more 
general conditions under which larger networks with con- 
tinuous dynamics are conform with their Boolean coun- 
terpart. 
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